module hope_recon_mod

  use spark_cubed_sphere
  use recon
  use params_mod
  use state_mod
  use diag_mod
  use static_mod

  implicit none

  private

  public hope_recon_init
  public find_trouble_cells
  public reconstruct
  public reconstruct_hor_ver
  public reconstruct_hor
  public reconstruct_ver
  public reconstruct_topo
  public reconstruct_topo_deriv
  public reconstruct_metric
  public hope_recon_final
  public recon_arrays
  public convolution
  public deconvolution
  public deconvolution_hor
  
  public rbs, rbe
  public rrs, rre
  public rts, rte
  public rls, rle
  public rqs, rqe
  public rds, rde
  public rus, rue

  ! There are few WENO types, so each has one instance.
  type(recon_type), allocatable, target :: recon_3d_obj(:)
  type(recon_type), allocatable, target :: recon_hor_obj(:)
  type(recon_type), allocatable, target :: recon_ver_obj(:)
  type(recon_type), allocatable, target :: recon_topo_obj(:)

  type recon_ptr_type
    type(recon_type), pointer :: ptr
  end type recon_ptr_type

  type recon_array_type
    type(recon_ptr_type), allocatable :: recon_3d (:,:,:)
    type(recon_ptr_type), allocatable :: recon_hor(:,:)
    type(recon_ptr_type), allocatable :: recon_ver(:)
    type(recon_ptr_type), allocatable :: recon_topo(:,:)
  end type recon_array_type

  ! Each cell on each domain has a pointer to the reconstruciton object instance.
  type(recon_array_type), allocatable :: recon_arrays(:)
  
  interface 
    subroutine reconstruct_3d_interface(state, diag)
      import state_type
      import diag_type
      type(state_type), intent(inout) :: state
      type(diag_type ), intent(inout) :: diag
    end subroutine reconstruct_3d_interface
  end interface

  procedure(reconstruct_3d_interface), pointer :: reconstruct => null()

  integer :: rbs, rbe
  integer :: rrs, rre
  integer :: rts, rte
  integer :: rls, rle
  integer :: rds, rde
  integer :: rus, rue
  integer :: rqs, rqe

contains

  subroutine hope_recon_init(mesh)

    type(cubed_sphere_mesh_type), intent(in) :: mesh

    integer i, j, k, ierr, iType
    integer idom, its, ite, jts, jte, kds, kde
    integer nd, n_recon_type_hor
    integer kc, irw
    real(r16)    xeq(nfq)
    real(r_kind) weq(nfq)
    
    integer, allocatable :: cell_mask(:,:,:)

    call log_notice('Initialize horizontal reconstruction functions.', pid=proc%id)
    
    ! Transfer coordinate from [-1, 1] to [-0.5, 0.5].
    call gaussian_legendre(nfq, xeq, weq)
    xeq = xeq / 2.0_16
    weq = weq / 2.0_16

    ! Set reconstruction target indices
    rbs = 1
    rbe = rbs + nfq - 1
    rrs = rbe + 1
    rre = rrs + nfq - 1
    rts = rre + 1
    rte = rts + nfq - 1
    rls = rte + 1
    rle = rls + nfq - 1
    rqs = rle + 1
    rqe = rqs + nvq - 1
    rds = rqe + 1
    rde = rds + nfq - 1
    rus = rde + 1
    rue = rus + nfq - 1
    
    ! Init hor recon structure for topography
    allocate(recon_topo_obj(0:0))
    recon_topo_obj(0)%init              => poly_square_init
    recon_topo_obj(0)%calc_recon_matrix => poly_square_calc_recon_matrix     
    recon_topo_obj(0)%reconstruct       => poly_square_reconstruct           
    recon_topo_obj(0)%reconstruct_deriv => poly_square_reconstruct_deriv
    
    ! Initialize each recon type with local coordinate origin at stencil center cell.
    call recon_topo_obj(0)%init(nd=2, sw=sw, dx=real(dx,16), dy=real(dy,16))
    
    ! Add reconstruction target points
    ! NOTE: The order matters.
    ! Bottom edge
    do k = 1, merge(neq, 1, hdim==3)
      do i = 1, neq
        call recon_topo_obj(0)%add_point(x=xeq(i), y=-0.5_16)
      end do
    end do
    ! Right edge
    do k = 1, merge(neq, 1, hdim==3)
      do i = 1, neq
        call recon_topo_obj(0)%add_point(x= 0.5_16, y=xeq(i))
      end do
    end do
    ! Top edge
    do k = 1, merge(neq, 1, hdim==3)
      do i = 1, neq
        call recon_topo_obj(0)%add_point(x=xeq(i), y= 0.5_16)
      end  do
    end do
    ! Left edge
    do k = 1, merge(neq, 1, hdim==3)
      do i = 1, neq
        call recon_topo_obj(0)%add_point(x=-0.5_16, y=xeq(i))
      end do
    end do
    ! Cell quadrature
    do k = 1, merge(neq, 1, hdim==3)
      do j = 1, neq
        do i = 1, neq
          call recon_topo_obj(0)%add_point(x=xeq(i), y=xeq(j))
        end do
      end do
    end do
    ! Down edge
    do j = 1, neq
      do i = 1, neq
        call recon_topo_obj(0)%add_point(x=xeq(i), y=xeq(j))
      end do
    end do
    ! Up Edge
    do j = 1, neq
      do i = 1, neq
        call recon_topo_obj(0)%add_point(x=xeq(i), y=xeq(j))
      end do
    end do
    
    call recon_topo_obj(0)%calc_recon_matrix(ierr)
    if (ierr /= 0) call log_error('Failed to calculate topo reconstruction matrix!', pid=proc%id)
    call recon_topo_obj(0)%release_unused_memory()
    
    call log_notice('Link topo reconstruction objects to each cell.', pid=proc%id)
    if(.not.allocated(recon_arrays))allocate(recon_arrays(mesh%ndom))
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte)
        allocate(recon_arrays(idom)%recon_topo(its:ite,jts:jte))
        do j = jts, jte
          do i = its, ite
            recon_arrays(idom)%recon_topo(i,j)%ptr => recon_topo_obj(0)
          end do
        end do
      end if
    end do

    ! Init recon structure for state
    if(hdim==2)then

      allocate(recon_hor_obj(0:0))
      
      select case(trim(recon_h_scheme))
      case ('WENO2D')
        nd = 2
        recon_hor_obj(0)%init                   => weno_square_init
        recon_hor_obj(0)%calc_recon_matrix      => weno_square_calc_recon_matrix     
        recon_hor_obj(0)%calc_ideal_coefs       => weno_square_calc_ideal_coefs      
        recon_hor_obj(0)%reconstruct            => weno_square_reconstruct           
        recon_hor_obj(0)%reconstruct_deriv      => weno_square_reconstruct_deriv
        recon_hor_obj(0)%trouble_cell_indicator => weno_square_trouble_cell_indicator
      case('WENO2D_rhombus')
        nd = 2
        recon_hor_obj(0)%init                   => weno_rhombus_init
        recon_hor_obj(0)%calc_recon_matrix      => weno_rhombus_calc_recon_matrix     
        recon_hor_obj(0)%calc_ideal_coefs       => weno_rhombus_calc_ideal_coefs      
        recon_hor_obj(0)%reconstruct            => weno_rhombus_reconstruct           
        recon_hor_obj(0)%reconstruct_deriv      => weno_rhombus_reconstruct_deriv
        recon_hor_obj(0)%trouble_cell_indicator => weno_rhombus_trouble_cell_indicator
      case ('Poly2D')
        nd = 2
        recon_hor_obj(0)%init                   => poly_square_init
        recon_hor_obj(0)%calc_recon_matrix      => poly_square_calc_recon_matrix     
        recon_hor_obj(0)%reconstruct            => poly_square_reconstruct           
        recon_hor_obj(0)%reconstruct_deriv      => poly_square_reconstruct_deriv
      case ('Poly2D_rhombus')
        nd = 2
        recon_hor_obj(0)%init                   => poly_rhombus_init
        recon_hor_obj(0)%calc_recon_matrix      => poly_rhombus_calc_recon_matrix     
        recon_hor_obj(0)%reconstruct            => poly_rhombus_reconstruct          
        recon_hor_obj(0)%reconstruct_deriv      => poly_rhombus_reconstruct_deriv
      end select
    
      ! Initialize each recon type with local coordinate origin at stencil center cell.
      call recon_hor_obj(0)%init(nd=nd, sw=sw, dx=real(dx,16), dy=real(dy,16))
      
      ! Add reconstruction target points
      ! NOTE: The order matters.
      ! Bottom edge
      do i = 1, neq
        call recon_hor_obj(0)%add_point(x=xeq(i), y=-0.5_16)
      end do
      ! Right edge
      do i = 1, neq
        call recon_hor_obj(0)%add_point(x= 0.5_16, y=xeq(i))
      end do
      ! Top edge
      do i = 1, neq
        call recon_hor_obj(0)%add_point(x=xeq(i), y= 0.5_16)
      end do
      ! Left edge
      do i = 1, neq
        call recon_hor_obj(0)%add_point(x=-0.5_16, y=xeq(i))
      end do
      ! Cell quadrature
      do j = 1, neq
        do i = 1, neq
          call recon_hor_obj(0)%add_point(x=xeq(i), y=xeq(j))
        end do
      end do
      
      if( trim(recon_h_scheme)=='WENO1D' .or. trim(recon_h_scheme)=='WENO2D' .or. trim(recon_h_scheme)=='WENO2D_rhombus' )then
        call recon_hor_obj(0)%calc_ideal_coefs(ierr)
        if (ierr /= 0) call log_error('Failed to calculate WENO ideal coefficients!', pid=proc%id)
      elseif(trim(recon_h_scheme)=='Poly1D' .or. trim(recon_h_scheme)=='Poly2D' .or. trim(recon_h_scheme)=='Poly2D_rhombus' )then
        call recon_hor_obj(0)%calc_recon_matrix(ierr)
        if (ierr /= 0) call log_error('Failed to calculate reconstruction matrix!', pid=proc%id)
      endif
      call recon_hor_obj(0)%release_unused_memory()
      
      call log_notice('Link reconstruction objects to each cell.', pid=proc%id)
      if(.not.allocated(recon_arrays))allocate(recon_arrays(mesh%ndom))
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%initialized) then
          call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte)
          allocate(recon_arrays(idom)%recon_hor(its:ite,jts:jte))
          do j = jts, jte
            do i = its, ite
              recon_arrays(idom)%recon_hor(i,j)%ptr => recon_hor_obj(0)
            end do
          end do
        end if
      end do
    elseif(hdim==3)then
      n_recon_type_hor = ( recon_h_order + 3 ) / 2
      allocate(recon_3d_obj(0:n_recon_type_hor-1))
      
      select case(trim(recon_h_scheme))
      case ('Poly3D')
        nd = 3
        do iType = 0, n_recon_type_hor - 1
          recon_3d_obj(iType)%init               => poly_cubic_init
          recon_3d_obj(iType)%calc_recon_matrix  => poly_cubic_calc_recon_matrix     
          recon_3d_obj(iType)%reconstruct        => poly_cubic_reconstruct          
          recon_3d_obj(iType)%reconstruct_deriv  => poly_cubic_reconstruct_deriv
        end do

        ! Initialize recon structure
        allocate( cell_mask(recon_h_order,recon_h_order,recon_h_order) )
        kc = ( recon_h_order + 1 ) / 2
        ! Center
        cell_mask = 1
        call recon_3d_obj(0)%init(nd=3, sw=sw, dx=real(dx,r16), dy=real(dy,r16), dz=real(deta,r16))
        ! Bottom
        cell_mask = 0
        cell_mask(:,:,kc:kc+2) = 1
        call recon_3d_obj(1)%init(nd=3, sw=sw, dx=real(dx,r16), dy=real(dy,r16), dz=real(deta,r16), mask=cell_mask)
        ! Top
        cell_mask = 0
        cell_mask(:,:,kc-2:kc) = 1
        call recon_3d_obj(2)%init(nd=3, sw=sw, dx=real(dx,r16), dy=real(dy,r16), dz=real(deta,r16), mask=cell_mask)
        ! Near boundary
        do irw = 1, rwv-1
          cell_mask = 0
          cell_mask(:,:,kc-irw:kc+irw) = 1
          call recon_3d_obj(2+irw)%init(nd=3, sw=sw, dx=real(dx,r16), dy=real(dy,r16), dz=real(deta,r16), mask=cell_mask)
        enddo
      
        ! Add reconstruction target points
        ! NOTE: The order matters.
        ! Bottom face
        do k = 1, neq
          do i = 1, neq
            do iType = 0, n_recon_type_hor - 1
              call recon_3d_obj(iType)%add_point(x=xeq(i), y=-0.5_16, z=xeq(k))
            end do
          end do
        end do
        ! Right face
        do k = 1, neq
          do j = 1, neq
            do iType = 0, n_recon_type_hor - 1
              call recon_3d_obj(iType)%add_point(x= 0.5_16, y=xeq(j), z=xeq(k))
            end do
          end do
        end do
        ! Top face
        do k = 1, neq
          do i = 1, neq
            do iType = 0, n_recon_type_hor - 1
              call recon_3d_obj(iType)%add_point(x=xeq(i), y= 0.5_16, z=xeq(k))
            end do
          enddo
        end do
        ! Left face
        do k = 1, neq
          do j = 1, neq
            do iType = 0, n_recon_type_hor - 1
              call recon_3d_obj(iType)%add_point(x=-0.5_16, y=xeq(j), z=xeq(k))
            end do
          end do
        end do
        ! Cell quadrature
        do k = 1, neq
          do j = 1, neq
            do i = 1, neq
              do iType = 0, n_recon_type_hor - 1
                call recon_3d_obj(iType)%add_point(x=xeq(i), y=xeq(j), z=xeq(k))
              end do
            end do
          end do
        end do
        ! Down face
        do j = 1, neq
          do i = 1, neq
            do iType = 0, n_recon_type_hor - 1
              call recon_3d_obj(iType)%add_point(x=xeq(i), y=xeq(j), z=-0.5_16)
            end do
          end do
        end do
        ! Up face
        do j = 1, neq
          do i = 1, neq
            do iType = 0, n_recon_type_hor - 1
              call recon_3d_obj(iType)%add_point(x=xeq(i), y=xeq(j), z=0.5_16)
            end do
          end do
        end do

        do iType = 0, n_recon_type_hor - 1
          call recon_3d_obj(iType)%calc_recon_matrix(ierr)
        enddo
        if (ierr /= 0) call log_error('Failed to calculate 3D reconstruction matrix!', pid=proc%id)
      
        call log_notice('Link 3D reconstruction objects to each cell.', pid=proc%id)
        if(.not.allocated(recon_arrays))allocate(recon_arrays(mesh%ndom))
        do idom = 1, mesh%ndom
          if (mesh%domain(idom)%initialized) then
            call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte, kds=kds, kde=kde)
            allocate(recon_arrays(idom)%recon_3d(its:ite,jts:jte,kds:kde))
            do j = jts, jte
              do i = its, ite
                k = kds
                recon_arrays(idom)%recon_3d(i,j,k)%ptr => recon_3d_obj(1)
                k = kde
                recon_arrays(idom)%recon_3d(i,j,k)%ptr => recon_3d_obj(2)
                do irw = 1, rwv-2
                  k = kds + irw
                  recon_arrays(idom)%recon_3d(i,j,k)%ptr => recon_3d_obj(2+irw)
                  k = kde - irw
                  recon_arrays(idom)%recon_3d(i,j,k)%ptr => recon_3d_obj(2+irw)
                enddo
              end do
            end do
          end if
        end do

      case default
        print*,'Only support recon_h_scheme="Poly3D" as reconstruction scheme when hdim==3'
        stop 'Check recon_h_scheme and hdim'
      end select

    endif
    
    if( model_type == 2 )then
      call log_notice('Calculate vertical reconstruction matrices.', pid=proc%id)
      
      if(vdim==1)then
        allocate(recon_ver_obj(0:3))
        
        iType = 0
        select case(trim(recon_v_scheme))
        case ('Poly1D')
          recon_ver_obj(iType)%init              => poly_1d_init
          recon_ver_obj(iType)%calc_recon_matrix => poly_1d_calc_recon_matrix     
          recon_ver_obj(iType)%reconstruct       => poly_1d_reconstruct           
          recon_ver_obj(iType)%reconstruct_deriv => poly_1d_reconstruct_deriv
          recon_ver_obj(iType)%get_jab           => poly_1d_get_jacobian
        case('WENO1D')                               
          recon_ver_obj(iType)%init              => weno_1d_init
          recon_ver_obj(iType)%calc_recon_matrix => weno_1d_calc_recon_matrix     
          recon_ver_obj(iType)%reconstruct       => weno_1d_reconstruct           
          recon_ver_obj(iType)%reconstruct_deriv => weno_1d_reconstruct_deriv
          !recon_ver_obj(iType)%get_jab           => weno_1d_get_jacobian ! WENO1D set get_jab in weno_1d_mod  
          recon_ver_obj(iType)%calc_ideal_coefs  => weno_1d_calc_ideal_coefs   
        case default
          stop 'Unknown recon_v_scheme'
        end select
        
        do iType = 1,3
          recon_ver_obj(iType)%init              => poly_1d_init
          recon_ver_obj(iType)%calc_recon_matrix => poly_1d_calc_recon_matrix     
          recon_ver_obj(iType)%reconstruct       => poly_1d_reconstruct           
          recon_ver_obj(iType)%reconstruct_deriv => poly_1d_reconstruct_deriv 
          recon_ver_obj(iType)%get_jab           => poly_1d_get_jacobian
        end do
        
        !if( recon_v_order /= 5 ) stop 'HOPE support only recon_v_order==5 for now'
        allocate( cell_mask(recon_v_order,1,1) )
        cell_mask(:,1,1) = [1,1,1,1,1]
        call recon_ver_obj(0)%init(nd=1, sw=swv, dx=real(deta,r16))
        cell_mask(:,1,1) = [0,0,1,1,1]
        call recon_ver_obj(1)%init(nd=1, sw=5, dx=real(deta,r16), mask=cell_mask)
        cell_mask(:,1,1) = [1,1,1,0,0]            
        call recon_ver_obj(2)%init(nd=1, sw=5, dx=real(deta,r16), mask=cell_mask)
        cell_mask(:,1,1) = [0,1,1,1,0]            
        call recon_ver_obj(3)%init(nd=1, sw=3, dx=real(deta,r16), mask=cell_mask(2:4,1:1,1:1))
        
        do iType = 0, 3
          ! Do not change the sequence of adding points, -0,5 must be the fist point, 0,5 must be the second one,
          ! or the Jacobian may not be calculated correctly.
          call recon_ver_obj(iType)%add_point(x=-0.5_16)
          call recon_ver_obj(iType)%add_point(x= 0.5_16)
        end do
        
        do iType = 0, 3
          call recon_ver_obj(iType)%calc_recon_matrix(ierr)
          if (ierr /= 0) call log_error('Failed to calculate vertical reconstruction matrix!', pid=proc%id)
        end do
        
        if( trim(recon_v_scheme)=='WENO1D' ) then
          call recon_ver_obj(0)%calc_ideal_coefs(ierr)
          if (ierr /= 0) call log_error('Failed to calculate vertical WENO ideal coefficients!', pid=proc%id)
        endif
        do idom = 1, mesh%ndom
          if (mesh%domain(idom)%initialized) then
            call mesh%get_params(idom=idom, kds=kds, kde=kde)
            allocate( recon_arrays(idom)%recon_ver(kds:kde) )
            do k = kds, kde
              if( k == kds )then ! Bottom boundary
                recon_arrays(idom)%recon_ver(k)%ptr => recon_ver_obj(1)
              elseif( k == kde )then ! Top boundary
                recon_arrays(idom)%recon_ver(k)%ptr => recon_ver_obj(2)
              elseif( ( k <= rwv .and. k > kds ) .or.      &
                      ( kde - k + 1 <= rwv .and. k < kde ) )then ! Near boundary
                recon_arrays(idom)%recon_ver(k)%ptr => recon_ver_obj(3)
              else ! Center region
                recon_arrays(idom)%recon_ver(k)%ptr => recon_ver_obj(0)
              end if
            end do
          end if
        end do
      endif
    endif
    
    ! Set reconstruction sheme
    if(model_type==0)then
      reconstruct => reconstruct_hor
    elseif(model_type==1)then
      reconstruct => reconstruct_hor
    elseif(model_type==2)then
      if(vdim==1)then
        reconstruct => reconstruct_hor_ver
      elseif(vdim==3)then
        reconstruct => reconstruct_3d
      else
        print*,'Not supported vdim, you need to set vdim=1 or 3'
        stop
      endif
    endif
    
  end subroutine hope_recon_init

  subroutine find_trouble_cells(state, diag)
    type(state_type), intent(inout) :: state
    type(diag_type ), intent(inout) :: diag

    integer idom, its, ite, jts, jte, kds, kde, iVar, i, j, k

    associate (mesh => state%mesh)
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%initialized) then
        call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte, kds=kds, kde=kde)
        associate (tci   => diag%domain(idom)%tci       , &
                   recon => recon_arrays(idom)%recon_hor, &
                   dx    => mesh%domain(idom)%dx        , &
                   q     => state%domain(idom)%q )
        !$OMP PARALLEL DO COLLAPSE(4)
        do iVar = 1, nVar
          do k = kds, kde
            do j = jts, jte
              do i = its, ite
                tci(i,j,k,iVar) = recon(i,j)%ptr%trouble_cell_indicator(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), dx)
              end do
            end do
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end if
    end do
    end associate

  end subroutine find_trouble_cells

  subroutine reconstruct_3d(state, diag)
    type(state_type), intent(inout) :: state
    type(diag_type ), intent(inout) :: diag

    integer pc, pes(6) ,pee(6), pqs, pqe
    integer idom, its, ite, jts, jte, kds, kde, iVar, ingb, ierr
    integer rwL, rwR, rwB, rwT, rwD, rwU
    integer i, j, k
    real(r_kind) fo(nfq*4+nvq), fov(2*nfqv)

    associate (mesh => state%mesh)
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%active) then
        associate (domain => mesh%domain(idom)           , &
                   Jab    => mesh%domain(idom)%Jab       , &
                   q      => state%domain(idom)%q        , &
                   qL     => diag%domain(idom)%qL        , &
                   qR     => diag%domain(idom)%qR        , &
                   qB     => diag%domain(idom)%qB        , &
                   qT     => diag%domain(idom)%qT        , &
                   qQ     => diag%domain(idom)%qQ        , &
                   qD     => diag%domain(idom)%qD        , &
                   qU     => diag%domain(idom)%qU        , &
                   tci    => diag%domain(idom)%tci       , &
                   recon  => recon_arrays(idom)%recon_3d   )
        call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte, kds=kds, kde=kde, pc= pc ,pes=pes, pee=pee, pqs=pqs, pqe=pqe)
        
        ! Reconstruct target point values
        !$OMP PARALLEL
        !$OMP DO PRIVATE(fo, fov, ierr, rwL, rwR, rwB, rwT, rwD, rwU) COLLAPSE(4)
        do iVar = 1, nVar
          do k = kds, kde
            do j = jts, jte
              do i = its, ite
                rwL = recon(i,j,k)%ptr%rwL
                rwR = recon(i,j,k)%ptr%rwR
                rwB = recon(i,j,k)%ptr%rwB
                rwT = recon(i,j,k)%ptr%rwT
                rwD = recon(i,j,k)%ptr%rwD
                rwU = recon(i,j,k)%ptr%rwU
                call recon(i,j,k)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                qL(:,i,j,k,iVar) = fo(rls:rle) * Jab(pes(1):pee(1),i,j,k)
                qR(:,i,j,k,iVar) = fo(rrs:rre) * Jab(pes(2):pee(2),i,j,k)
                qB(:,i,j,k,iVar) = fo(rbs:rbe) * Jab(pes(3):pee(3),i,j,k)
                qT(:,i,j,k,iVar) = fo(rts:rte) * Jab(pes(4):pee(4),i,j,k)
                qD(:,i,j,k,iVar) = fo(rds:rde) * Jab(pes(5):pee(5),i,j,k)
                qU(:,i,j,k,iVar) = fo(rus:rue) * Jab(pes(6):pee(6),i,j,k)
                qQ(:,i,j,k,iVar) = fo(rqs:rqe) * Jab(pqs   :pqe   ,i,j,k)
              end do
            end do
          end do
        end do
        !$OMP END DO
        ! Fill vertical boundary
        !$OMP DO COLLAPSE(3)
        do iVar = 1, nVar
          do j = jts, jte
            do i = its, ite
              qU(:,i,j,kds-1,iVar) = qD(:,i,j,kds,iVar)
              qD(:,i,j,kde+1,iVar) = qU(:,i,j,kde,iVar)
            end do
          end do
        end do
        !$OMP END DO
        !$OMP END PARALLEL
        end associate
        
        ! Neighbor domains
        do ingb = 1, 4
          associate (ngb => mesh%domain(idom)%ngb(ingb))
          if (ngb%ptr%initialized .and. .not. ngb%ptr%active) then
            call mesh%get_params(idom=ngb%ptr%id, ids=its, ide=ite, jds=jts, jde=jte, kds=kds, kde=kde, pc= pc ,pes=pes, pee=pee, pqs=pqs, pqe=pqe)
            select case (ngb%orient)
            case (bottom)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qB    => diag%domain(ngb%ptr%id)%qB    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_3d)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j,k)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qB(:,i,j,k,iVar) = fo(rbs:rbe) * Jab(pes(3):pee(3),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (right)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qR    => diag%domain(ngb%ptr%id)%qR    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_3d)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j,k)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qR(:,i,j,k,iVar) = fo(rrs:rre) * Jab(pes(2):pee(2),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (top)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qT    => diag%domain(ngb%ptr%id)%qT    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_3d)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j,k)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qT(:,i,j,k,iVar) = fo(rts:rte) * Jab(pes(4):pee(4),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (left)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qL    => diag%domain(ngb%ptr%id)%qL    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_3d)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j,k)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qL(:,i,j,k,iVar) = fo(rls:rle) * Jab(pes(1):pee(1),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            end select
          end if
          end associate
        end do
      end if
    end do
    end associate

  end subroutine reconstruct_3d

  subroutine reconstruct_hor_ver(state, diag)
    type(state_type), intent(inout) :: state
    type(diag_type ), intent(inout) :: diag

    integer pc, pes(6) ,pee(6), pqs, pqe
    integer idom, its, ite, jts, jte, kds, kde, iVar, ingb, ierr
    integer rwL, rwR
    integer i, j, k
    real(r_kind) fo(nfq*4+nvq), fov(2*nfqv)

    associate (mesh => state%mesh)
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%active) then
        associate (domain => mesh%domain(idom)           , &
                   Jab    => mesh%domain(idom)%Jab       , &
                   q      => state%domain(idom)%q        , &
                   qL     => diag%domain(idom)%qL        , &
                   qR     => diag%domain(idom)%qR        , &
                   qB     => diag%domain(idom)%qB        , &
                   qT     => diag%domain(idom)%qT        , &
                   qQ     => diag%domain(idom)%qQ        , &
                   qD     => diag%domain(idom)%qD        , &
                   qU     => diag%domain(idom)%qU        , &
                   tci    => diag%domain(idom)%tci       , &
                   recon_h=> recon_arrays(idom)%recon_hor, &
                   recon_v=> recon_arrays(idom)%recon_ver )
        call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte, kds=kds, kde=kde, pc= pc ,pes=pes, pee=pee, pqs=pqs, pqe=pqe)
        
        ! Reconstruct target point values
        !$OMP PARALLEL
        !$OMP DO PRIVATE(fo, fov, ierr, rwL, rwR) COLLAPSE(4)
        do iVar = 1, nVar
          do k = kds, kde
            do j = jts, jte
              do i = its, ite
                call recon_h(i,j)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                qB(:,i,j,k,iVar) = fo(rbs:rbe) * Jab(pes(3):pee(3),i,j,k)
                qR(:,i,j,k,iVar) = fo(rrs:rre) * Jab(pes(2):pee(2),i,j,k)
                qT(:,i,j,k,iVar) = fo(rts:rte) * Jab(pes(4):pee(4),i,j,k)
                qL(:,i,j,k,iVar) = fo(rls:rle) * Jab(pes(1):pee(1),i,j,k)
                qQ(:,i,j,k,iVar) = fo(rqs:rqe) * Jab(pqs   :pqe   ,i,j,k)
                
                rwL = recon_v(k)%ptr%rwL
                rwR = recon_v(k)%ptr%rwR
                call recon_v(k)%ptr%reconstruct(q(i:i,j:j,k-rwL:k+rwR,iVar), fov)
                qD(:,i,j,k,iVar) = fov(1     :  nfqv) * Jab(pes(5):pee(5),i,j,k)
                qU(:,i,j,k,iVar) = fov(nfqv+1:2*nfqv) * Jab(pes(6):pee(6),i,j,k)
              end do
            end do
          end do
        end do
        !$OMP END DO
        ! Fill vertical boundary
        !$OMP DO COLLAPSE(3)
        do iVar = 1, nVar
          do j = jts, jte
            do i = its, ite
              qU(:,i,j,kds-1,iVar) = qD(:,i,j,kds,iVar)
              qD(:,i,j,kde+1,iVar) = qU(:,i,j,kde,iVar)
            end do
          end do
        end do
        !$OMP END DO
        !$OMP END PARALLEL
        end associate
        
        ! Neighbor domains
        do ingb = 1, 4
          associate (ngb => mesh%domain(idom)%ngb(ingb))
          if (ngb%ptr%initialized .and. .not. ngb%ptr%active) then
            call mesh%get_params(idom=ngb%ptr%id, ids=its, ide=ite, jds=jts, jde=jte, kds=kds, kde=kde, pc= pc ,pes=pes, pee=pee, pqs=pqs, pqe=pqe)
            select case (ngb%orient)
            case (bottom)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qB    => diag%domain(ngb%ptr%id)%qB    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qB(:,i,j,k,iVar) = fo(rbs:rbe) * Jab(pes(3):pee(3),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (right)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qR    => diag%domain(ngb%ptr%id)%qR    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qR(:,i,j,k,iVar) = fo(rrs:rre) * Jab(pes(2):pee(2),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (top)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qT    => diag%domain(ngb%ptr%id)%qT    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qT(:,i,j,k,iVar) = fo(rts:rte) * Jab(pes(4):pee(4),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (left)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qL    => diag%domain(ngb%ptr%id)%qL    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qL(:,i,j,k,iVar) = fo(rls:rle) * Jab(pes(1):pee(1),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            end select
          end if
          end associate
        end do
      end if
    end do
    end associate

  end subroutine reconstruct_hor_ver
  
  subroutine reconstruct_hor(state, diag)
    type(state_type), intent(inout) :: state
    type(diag_type ), intent(inout) :: diag

    integer pc, pes(6), pee(6), pqs, pqe, its, ite, jts, jte, kds, kde, iVar, ingb, ierr
    integer rwL, rwR
    integer idom, i, j, k
    real(r_kind) fo(nfq*4+nvq), fov(2*nfqv)

    associate (mesh => state%mesh)
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%active) then
        associate (domain => mesh%domain(idom)       , &
                   Jab    => mesh%domain(idom)%Jab   , &
                   q      => state%domain(idom)%q    , &
                   qL     => diag%domain(idom)%qL    , &
                   qR     => diag%domain(idom)%qR    , &
                   qB     => diag%domain(idom)%qB    , &
                   qT     => diag%domain(idom)%qT    , &
                   qQ     => diag%domain(idom)%qQ    , &
                   qD     => diag%domain(idom)%qD    , &
                   qU     => diag%domain(idom)%qU    , &
                   tci    => diag%domain(idom)%tci   , &
                   recon  => recon_arrays(idom)%recon_hor  )
        call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte, kds=kds, kde=kde, pc=pc, pes=pes, pee=pee, pqs=pqs, pqe=pqe)
        
        ! Reconstruct target point values
        !$OMP PARALLEL
        !$OMP DO PRIVATE(fo, ierr) COLLAPSE(4)
        do iVar = 1, nVar
          do k = kds, kde
            do j = jts, jte
              do i = its, ite
                call recon(i,j)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                qB(:,i,j,k,iVar) = fo(rbs:rbe) * Jab(pes(3):pee(3),i,j,k)
                qR(:,i,j,k,iVar) = fo(rrs:rre) * Jab(pes(2):pee(2),i,j,k)
                qT(:,i,j,k,iVar) = fo(rts:rte) * Jab(pes(4):pee(4),i,j,k)
                qL(:,i,j,k,iVar) = fo(rls:rle) * Jab(pes(1):pee(1),i,j,k)
                qQ(:,i,j,k,iVar) = fo(rqs:rqe) * Jab(pqs   :pqe   ,i,j,k)
              end do
            end do
          end do
        end do
        !$OMP END DO
        !$OMP END PARALLEL
        end associate
        
        ! Neighbor domains
        do ingb = 1, 4
          associate (ngb => mesh%domain(idom)%ngb(ingb))
          if (ngb%ptr%initialized .and. .not. ngb%ptr%active) then
            call mesh%get_params(idom=ngb%ptr%id, ids=its, ide=ite, jds=jts, jde=jte, kds=kds, kde=kde, pc=pc, pes=pes, pee=pee)
            select case (ngb%orient)
            case (bottom)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qB    => diag%domain(ngb%ptr%id)%qB    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qB(:,i,j,k,iVar) = fo(rbs:rbe) * Jab(pes(3):pee(3),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (right)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qR    => diag%domain(ngb%ptr%id)%qR    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qR(:,i,j,k,iVar) = fo(rrs:rre) * Jab(pes(2):pee(2),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (top)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qT    => diag%domain(ngb%ptr%id)%qT    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qT(:,i,j,k,iVar) = fo(rts:rte) * Jab(pes(4):pee(4),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (left)
              associate (Jab   => mesh%domain(idom)%Jab         , &
                         q     => state%domain(ngb%ptr%id)%q    , &
                         qL    => diag%domain(ngb%ptr%id)%qL    , &
                         tci   => diag%domain(ngb%ptr%id)%tci   , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      call recon(i,j)%ptr%reconstruct(q(i-rw:i+rw,j-rw:j+rw,k:k,iVar), fo, TCI=TCI(i,j,k,iVar))
                      qL(:,i,j,k,iVar) = fo(rls:rle) * Jab(pes(1):pee(1),i,j,k)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            end select
          end if
          end associate
        end do
      end if
    end do
    end associate

  end subroutine reconstruct_hor
  
  subroutine deconvolution_hor(state, diag)
    type(state_type), intent(inout) :: state
    type(diag_type ), intent(inout) :: diag

    integer idom, its, ite, jts, jte, kds, kde, iVar, ingb, pes(6), pee(6)
    integer rwL, rwR
    integer p, i, j, k
    real(r_kind) :: fxc(-1:nz+2), fyc(-1:nz+2) ! Use for convolution
    real(r_kind) :: fo(nfq*4+nvq), fov(2*nfqv)

    associate (mesh => state%mesh)
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%active) then
        associate (domain => mesh%domain(idom)   , &
                   qL     => diag%domain(idom)%qL, &
                   qR     => diag%domain(idom)%qR, &
                   qB     => diag%domain(idom)%qB, &
                   qT     => diag%domain(idom)%qT)
        call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte, kds=kds, kde=kde, pes=pes, pee=pee)
        
        ! Reconstruct target point values
        !$OMP PARALLEL
        !$OMP DO COLLAPSE(3) PRIVATE(p,k,fxc)
        do iVar = 1, nVar
          do j = jts, jte
            do i = its, ite
              do p = pes(1), pee(1)
                fxc(1:nz) = qL(p,i,j,kds:kde,iVar)
                do k = kds, kde
                  qL(p,i,j,k,iVar) = deconvolution( fxc(k-2:k+2), k, kds, kde )
                enddo
              enddo
              do p = pes(2), pee(2)
                fxc(1:nz) = qR(p,i,j,kds:kde,iVar)
                do k = kds, kde
                  qR(p,i,j,k,iVar) = deconvolution( fxc(k-2:k+2), k, kds, kde )
                enddo
              enddo
            enddo
          enddo
        enddo
        !$OMP END DO NOWAIT
        
        !$OMP DO COLLAPSE(3) PRIVATE(p,k,fyc)
        do iVar = 1, nVar
          do j = jts, jte
            do i = its, ite
              do p = pes(3), pee(3)
                fyc(1:nz) = qB(p,i,j,kds:kde,iVar)
                do k = kds, kde
                  qB(p,i,j,k,iVar) = deconvolution( fyc(k-2:k+2), k, kds, kde )
                enddo
              enddo
              do p = pes(4), pee(4)
                fyc(1:nz) = qT(p,i,j,kds:kde,iVar)
                do k = kds, kde
                  qT(p,i,j,k,iVar) = deconvolution( fyc(k-2:k+2), k, kds, kde )
                enddo
              enddo
            enddo
          enddo
        enddo
        !$OMP END DO
        !$OMP END PARALLEL
        end associate
        
        ! Neighbor domains
        do ingb = 1, 4
          associate (ngb => mesh%domain(idom)%ngb(ingb))
          if (ngb%ptr%initialized .and. .not. ngb%ptr%active) then
            call mesh%get_params(idom=ngb%ptr%id, ids=its, ide=ite, jds=jts, jde=jte, kds=kds, kde=kde)
            select case (ngb%orient)
            case (bottom)
              associate (qB => diag%domain(ngb%ptr%id)%qB)
              !$OMP PARALLEL DO PRIVATE(fo, p, i, j, k) COLLAPSE(4)
              do iVar = 1, nVar
                do j = jts, jte
                  do i = its, ite
                    do p = pes(3), pee(3)
                      fyc(1:nz) = qB(p,i,j,kds:kde,iVar)
                      do k = kds, kde
                        qB(p,i,j,k,iVar) = deconvolution( fyc(k-2:k+2), k, kds, kde )
                      enddo
                    enddo
                  enddo
                enddo
              end do
              !$OMP END PARALLEL DO
              end associate
            case (right)
              associate (qR => diag%domain(ngb%ptr%id)%qR)
              !$OMP PARALLEL DO PRIVATE(fo, p, i, j, k) COLLAPSE(4)
              do iVar = 1, nVar
                do j = jts, jte
                  do i = its, ite
                    do p = pes(2), pee(2)
                      fxc(1:nz) = qR(p,i,j,kds:kde,iVar)
                      do k = kds, kde
                        qR(p,i,j,k,iVar) = deconvolution( fxc(k-2:k+2), k, kds, kde )
                      enddo
                    enddo
                  enddo
                enddo
              end do
              !$OMP END PARALLEL DO
              end associate
            case (top)
              associate (qT => diag%domain(ngb%ptr%id)%qT)
              !$OMP PARALLEL DO PRIVATE(fo, p, i, j, k) COLLAPSE(4)
              do iVar = 1, nVar
                do j = jts, jte
                  do i = its, ite
                    do p = pes(4), pee(4)
                      fyc(1:nz) = qT(p,i,j,kds:kde,iVar)
                      do k = kds, kde
                        qT(p,i,j,k,iVar) = deconvolution( fyc(k-2:k+2), k, kds, kde )
                      enddo
                    enddo
                  enddo
                enddo
              end do
              !$OMP END PARALLEL DO
              end associate
            case (left)
              associate (qL => diag%domain(ngb%ptr%id)%qL)
              !$OMP PARALLEL DO PRIVATE(fo, p, i, j, k) COLLAPSE(4)
              do iVar = 1, nVar
                do j = jts, jte
                  do i = its, ite
                    do p = pes(1), pee(1)
                      fxc(1:nz) = qL(p,i,j,kds:kde,iVar)
                      do k = kds, kde
                        qL(p,i,j,k,iVar) = deconvolution( fxc(k-2:k+2), k, kds, kde )
                      enddo
                    enddo
                  enddo
                enddo
              end do
              !$OMP END PARALLEL DO
              end associate
            end select
          end if
          end associate
        end do
      end if
    end do
    end associate

  end subroutine deconvolution_hor

  subroutine reconstruct_ver(state, diag, idom, i, j)
    type(state_type), intent(inout) :: state
    type(diag_type ), intent(inout) :: diag
    integer,          intent(in   ) :: idom,i,j
    
    integer pc, pes(6), pee(6), kds, kde
    integer iVar, k, nc
    integer rwL, rwR
    real(r_kind) fov(2)
    real(r_kind) qc(swv,1,1)

    associate (mesh    => state%mesh,                 &
               Jab     => mesh%domain(idom)%Jab,      &
               q       => state%domain(idom)%q,       &
               qD      => diag%domain(idom)%qD,       &
               qU      => diag%domain(idom)%qU,       &
               recon_v => recon_arrays(idom)%recon_ver)
    call mesh%get_params(idom=idom, pc=pc,pes=pes, pee=pee, kds=kds, kde=kde)
    do iVar = 1, nVar
      do k = kds, kde
        rwL = recon_v(k)%ptr%rwL
        rwR = recon_v(k)%ptr%rwR
        nc  = rwR + rwL + 1
        
        !qc(1:nc,1,1) = q(i,j,k-rwL:k+rwR,iVar)
        !call recon_v(k)%ptr%reconstruct(qc(1:nc,:,:), fov)
        !qD(:,i,j,k,iVar) = fov(1)
        !qU(:,i,j,k,iVar) = fov(2)
        
        qc(1:nc,1,1) = q(i,j,k-rwL:k+rwR,iVar)
        call recon_v(k)%ptr%reconstruct(qc(1:nc,:,:), fov)
        qD(:,i,j,k,iVar) = fov(1) * Jab(pes(5):pee(5),i,j,k)
        qU(:,i,j,k,iVar) = fov(2) * Jab(pes(6):pee(6),i,j,k)
      end do
      ! Fill boundary
      qU(:,i,j,kds-1,iVar) = qD(:,i,j,kds,iVar)
      qD(:,i,j,kde+1,iVar) = qU(:,i,j,kde,iVar)
    end do
    end associate
  end subroutine reconstruct_ver
  
  subroutine reconstruct_topo(static)
    type(static_type), intent(inout) :: static

    integer idom, its, ite, jts, jte, kds, kde, iVar, ingb, ierr
    integer i, j, k
    real(r_kind) fo(nfq*4+nvq+2*nfqv)
    real(r_kind) q(-rw:rw,-rw:rw,1)
    
    associate( mesh => static%mesh )
      ! Topo interpolate to points
      do idom = 1, mesh%ndom
        if (.not. mesh%domain(idom)%initialized) cycle
        associate (domain => mesh%domain(idom),      &
                   zs     => static%domain(idom)%zs, &
                   zsL    => static%domain(idom)%zsL,&
                   zsR    => static%domain(idom)%zsR,&
                   zsB    => static%domain(idom)%zsB,&
                   zsT    => static%domain(idom)%zsT,&
                   zsD    => static%domain(idom)%zsD,&
                   zsU    => static%domain(idom)%zsU,&
                   zsQ    => static%domain(idom)%zsQ,&
                   recon  => recon_arrays(idom)%recon_topo )
        call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte)
        ! Reconstruct target point values
        !$OMP PARALLEL DO PRIVATE(q, fo) COLLAPSE(2)
        do j = jts, jte
          do i = its, ite
            q(:,:,1) = zs(i-rw:i+rw,j-rw:j+rw)
            call recon(i,j)%ptr%reconstruct(q(:,:,1:1), fo, TCI=.false.)
            zsB(:,i,j) = fo(rbs:rbe)
            zsR(:,i,j) = fo(rrs:rre)
            zsT(:,i,j) = fo(rts:rte)
            zsL(:,i,j) = fo(rls:rle)
            zsD(:,i,j) = fo(rds:rde)
            zsU(:,i,j) = fo(rus:rue)
            zsQ(:,i,j) = fo(rqs:rqe)
          end do
        end do
        !$OMP END PARALLEL DO
        
        ! Neighbor domains
        do ingb = 1, 4
          associate (ngb => mesh%domain(idom)%ngb(ingb))
          if (ngb%ptr%initialized .and. .not. ngb%ptr%active) then
            call mesh%get_params(idom=ngb%ptr%id, ids=its, ide=ite, jds=jts, jde=jte, kds=kds, kde=kde)
            select case (ngb%orient)
            case (bottom)
              associate (zs    => static%domain(idom)%zs,  &
                         zsB   => static%domain(idom)%zsB, &
                         recon => recon_arrays(ngb%ptr%id)%recon_topo)
              !$OMP PARALLEL DO PRIVATE(q, fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      q(:,:,1) = zs(i-rw:i+rw,j-rw:j+rw)
                      call recon(i,j)%ptr%reconstruct(q(:,:,1:1), fo, TCI=.false.)
                      zsB(:,i,j) = fo(rbs:rbe)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (right)
              associate (zs    => static%domain(idom)%zs,  &
                         zsR   => static%domain(idom)%zsR, &
                         recon => recon_arrays(ngb%ptr%id)%recon_topo)
              !$OMP PARALLEL DO PRIVATE(q, fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      q(:,:,1) = zs(i-rw:i+rw,j-rw:j+rw)
                      call recon(i,j)%ptr%reconstruct(q(:,:,1:1), fo, TCI=.false.)
                      zsR(:,i,j) = fo(rrs:rre)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (top)
              associate (zs    => static%domain(idom)%zs,  &
                         zsT   => static%domain(idom)%zsT, &
                         recon => recon_arrays(ngb%ptr%id)%recon_topo)
              !$OMP PARALLEL DO PRIVATE(q, fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      q(:,:,1) = zs(i-rw:i+rw,j-rw:j+rw)
                      call recon(i,j)%ptr%reconstruct(q(:,:,1:1), fo, TCI=.false.)
                      zsT(:,i,j) = fo(rts:rte)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (left)
              associate (zs    => static%domain(idom)%zs,  &
                         zsL   => static%domain(idom)%zsL, &
                         recon => recon_arrays(ngb%ptr%id)%recon_topo)
              !$OMP PARALLEL DO PRIVATE(q, fo, ierr) COLLAPSE(4)
              do iVar = 1, nVar
                do k = kds, kde
                  do j = jts, jte
                    do i = its, ite
                      q(:,:,1) = zs(i-rw:i+rw,j-rw:j+rw)
                      call recon(i,j)%ptr%reconstruct(q(:,:,1:1), fo, TCI=.false.)
                      zsL(:,i,j) = fo(rls:rle)
                    end do
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            end select
          end if
          end associate
        end do
        end associate
      end do
    end associate
  end subroutine reconstruct_topo
  
  subroutine reconstruct_topo_deriv(static)
    type(static_type), intent(inout) :: static

    integer idom, its, ite, jts, jte, kds, kde, iVar, ingb, ierr, i, j, k
    real(r_kind) fx(nfq*4+nvq+2*nfqv), fy(nfq*4+nvq+2*nfqv)
    real(r_kind) q(-rw:rw,-rw:+rw,1)
    
    associate (mesh => static%mesh)
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%active) then
          associate (domain => mesh%domain(idom)         , &
                     zs     => static%domain(idom)%zs    , &
                     dzsdxL => static%domain(idom)%dzsdxL, &
                     dzsdxR => static%domain(idom)%dzsdxR, &
                     dzsdxB => static%domain(idom)%dzsdxB, &
                     dzsdxT => static%domain(idom)%dzsdxT, &
                     dzsdxD => static%domain(idom)%dzsdxD, &
                     dzsdxU => static%domain(idom)%dzsdxU, &
                     dzsdxQ => static%domain(idom)%dzsdxQ, &
                     dzsdyL => static%domain(idom)%dzsdyL, &
                     dzsdyR => static%domain(idom)%dzsdyR, &
                     dzsdyB => static%domain(idom)%dzsdyB, &
                     dzsdyT => static%domain(idom)%dzsdyT, &
                     dzsdyD => static%domain(idom)%dzsdyD, &
                     dzsdyU => static%domain(idom)%dzsdyU, &
                     dzsdyQ => static%domain(idom)%dzsdyQ, &
                     recon  => recon_arrays(idom)%recon_topo)
          call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte, kds=kds, kde=kde)
          
          ! Reconstruct target point values
          !$OMP PARALLEL DO PRIVATE(q, fx, fy, ierr) COLLAPSE(2)
          do j = jts, jte
            do i = its, ite
              q(:,:,1) = zs(i-rw:i+rw,j-rw:j+rw)
              call recon(i,j)%ptr%reconstruct_deriv(q(:,:,1:1), fx, fy)
              dzsdxB(:,i,j) = fx(rbs:rbe)
              dzsdxR(:,i,j) = fx(rrs:rre)
              dzsdxT(:,i,j) = fx(rts:rte)
              dzsdxL(:,i,j) = fx(rls:rle)
              dzsdxD(:,i,j) = fx(rds:rde)
              dzsdxU(:,i,j) = fx(rus:rue)
              dzsdxQ(:,i,j) = fx(rqs:rqe)
              dzsdyB(:,i,j) = fy(rbs:rbe)
              dzsdyR(:,i,j) = fy(rrs:rre)
              dzsdyT(:,i,j) = fy(rts:rte)
              dzsdyL(:,i,j) = fy(rls:rle)
              dzsdyD(:,i,j) = fy(rds:rde)
              dzsdyU(:,i,j) = fy(rus:rue)
              dzsdyQ(:,i,j) = fy(rqs:rqe)
            end do
          end do
          !$OMP END PARALLEL DO
          end associate
          
          ! Neighbor domains
          do ingb = 1, 4
            associate (ngb => mesh%domain(idom)%ngb(ingb))
            if (ngb%ptr%initialized .and. .not. ngb%ptr%active) then
              call mesh%get_params(idom=ngb%ptr%id, ids=its, ide=ite, jds=jts, jde=jte, kds=kds, kde=kde)
              select case (ngb%orient)
              case (bottom)
                associate (zs     => static%domain(ngb%ptr%id)%zs    , &
                           dzsdxB => static%domain(ngb%ptr%id)%dzsdxB, &
                           dzsdyB => static%domain(ngb%ptr%id)%dzsdyB, &
                           recon  => recon_arrays(ngb%ptr%id)%recon_topo)
                !$OMP PARALLEL DO PRIVATE(q, fx, fy, ierr) COLLAPSE(2)
                do j = jts, jte
                  do i = its, ite
                    q(:,:,1) = zs(i-rw:i+rw,j-rw:j+rw)
                    call recon(i,j)%ptr%reconstruct_deriv(q(:,:,1:1), fx, fy)
                    dzsdxB(:,i,j) = fx(rbs:rbe)
                    dzsdyB(:,i,j) = fy(rbs:rbe)
                  end do
                end do
                !$OMP END PARALLEL DO
                end associate
              case (right)
                associate (zs     => static%domain(ngb%ptr%id)%zs    , &
                           dzsdxR => static%domain(ngb%ptr%id)%dzsdxR, &
                           dzsdyR => static%domain(ngb%ptr%id)%dzsdyR, &
                           recon  => recon_arrays(ngb%ptr%id)%recon_topo   )
                !$OMP PARALLEL DO PRIVATE(q, fx, fy, ierr) COLLAPSE(2)
                do j = jts, jte
                  do i = its, ite
                    q(:,:,1) = zs(i-rw:i+rw,j-rw:j+rw)
                    call recon(i,j)%ptr%reconstruct_deriv(q(:,:,1:1), fx, fy)
                    dzsdxR(:,i,j) = fx(rrs:rre)
                    dzsdyR(:,i,j) = fy(rrs:rre)
                  end do
                end do
                !$OMP END PARALLEL DO
                end associate
              case (top)
                associate (zs     => static%domain(ngb%ptr%id)%zs     , &
                           dzsdxT => static%domain(ngb%ptr%id)%dzsdxT , &
                           dzsdyT => static%domain(ngb%ptr%id)%dzsdyT , &
                           recon  => recon_arrays(ngb%ptr%id)%recon_topo)
                !$OMP PARALLEL DO PRIVATE(q, fx, fy, ierr) COLLAPSE(2)
                do j = jts, jte
                  do i = its, ite
                    q(:,:,1) = zs(i-rw:i+rw,j-rw:j+rw)
                    call recon(i,j)%ptr%reconstruct_deriv(q(:,:,1:1), fx, fy)
                    dzsdxT(:,i,j) = fx(rts:rte)
                    dzsdyT(:,i,j) = fy(rts:rte)
                  end do
                end do
                !$OMP END PARALLEL DO
                end associate
              case (left)
                associate (zs     => static%domain(ngb%ptr%id)%zs     , &
                           dzsdxL => static%domain(ngb%ptr%id)%dzsdxL , &
                           dzsdyL => static%domain(ngb%ptr%id)%dzsdyL , &
                           recon => recon_arrays(ngb%ptr%id)%recon_topo)
                !$OMP PARALLEL DO PRIVATE(q, fx, fy, ierr) COLLAPSE(2)
                do j = jts, jte
                  do i = its, ite
                    q(:,:,1) = zs(i-rw:i+rw,j-rw:j+rw)
                    call recon(i,j)%ptr%reconstruct_deriv(q(:,:,1:1), fx, fy)
                    dzsdxL(:,i,j) = fx(rls:rle)
                    dzsdyL(:,i,j) = fy(rls:rle)
                  end do
                end do
                !$OMP END PARALLEL DO
                end associate
              end select
            end if
            end associate
          end do
        end if
      end do
    end associate
  end subroutine reconstruct_topo_deriv

  subroutine reconstruct_metric(static)

    type(static_type), intent(inout) :: static

    integer idom, pc, pes(6), pee(6), pqs, pqe, its, ite, jts, jte, kds, kde, iVar, ingb, ierr
    integer rwL, rwR
    integer i, j, k, icol, irow
    real(r_kind) fo(nfq*4+nvq), fov(2*nfqv)

    associate (mesh => static%mesh)
    do idom = 1, mesh%ndom
      if (mesh%domain(idom)%active) then
        associate (domain => mesh%domain(idom)           , &
                   q      => mesh%domain(idom)%Jab       , &
                   recon_h=> recon_arrays(idom)%recon_hor, &
                   recon_v=> recon_arrays(idom)%recon_ver)
        call mesh%get_params(idom=idom, pc=pc, pes=pes, pee=pee, pqs=pqs, pqe=pqe,&
                             its=its, ite=ite, jts=jts, jte=jte, kds=kds, kde=kde)
        
        ! Reconstruct target point values
        !$OMP PARALLEL
        !$OMP DO PRIVATE(fo, fov, ierr, rwL, rwR) COLLAPSE(3)
        do k = kds, kde
          do j = jts, jte
            do i = its, ite
              call recon_h(i,j)%ptr%reconstruct(q(pc,i-rw:i+rw,j-rw:j+rw,k:k), fo, TCI=.false.)
              q(pes(1):pee(1),i,j,k) = fo(rls:rle)
              q(pes(2):pee(2),i,j,k) = fo(rrs:rre)
              q(pes(3):pee(3),i,j,k) = fo(rbs:rbe)
              q(pes(4):pee(4),i,j,k) = fo(rts:rte)
              q(pqs   :pqe   ,i,j,k) = fo(rqs:rqe)
              
              rwL = recon_v(k)%ptr%rwL
              rwR = recon_v(k)%ptr%rwR
              call recon_v(k)%ptr%reconstruct(q(i:i,j:j,k-rwL:k+rwR,iVar), fov)
              q(pes(5):pee(5),i,j,k) = fov(1:nfqv)
              q(pes(6):pee(6),i,j,k) = fov(nfqv+1:2*nfqv)
            end do
          end do
        end do
        !$OMP END DO
        ! Fill vertical boundary
        !$OMP DO COLLAPSE(2)
        do j = jts, jte
          do i = its, ite
            q(pes(6):pee(6),i,j,kds-1) = q(pes(5):pee(5),i,j,kds)
            q(pes(5):pee(5),i,j,kde+1) = q(pes(6):pee(6),i,j,kde)
          end do
        end do
        !$OMP END DO
        !$OMP END PARALLEL
        end associate
        
        ! Neighbor domains
        do ingb = 1, 4
          associate (ngb => mesh%domain(idom)%ngb(ingb))
          if (ngb%ptr%initialized .and. .not. ngb%ptr%active) then
            call mesh%get_params(idom=ngb%ptr%id, pc=pc, pes=pes, pee=pee, &
                                 ids=its, ide=ite, jds=jts, jde=jte, kds=kds, kde=kde)
            select case (ngb%orient)
            case (bottom)
              associate (q     => mesh%domain(ngb%ptr%id)%Jab    , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(3)
              do k = kds, kde
                do j = jts, jte
                  do i = its, ite
                    call recon(i,j)%ptr%reconstruct(q(pc,i-rw:i+rw,j-rw:j+rw,k:k), fo, TCI=.false.)
                    q(pes(3):pee(3),i,j,k) = fo(rbs:rbe)
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (right)
              associate (q     => mesh%domain(ngb%ptr%id)%Jab  , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(3)
              do k = kds, kde
                do j = jts, jte
                  do i = its, ite
                    call recon(i,j)%ptr%reconstruct(q(pc,i-rw:i+rw,j-rw:j+rw,k:k), fo, TCI=.false.)
                    q(pes(2):pee(2),i,j,k) = fo(rrs:rre)
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (top)
              associate (q     => mesh%domain(ngb%ptr%id)%Jab  , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(3)
              do k = kds, kde
                do j = jts, jte
                  do i = its, ite
                    call recon(i,j)%ptr%reconstruct(q(pc,i-rw:i+rw,j-rw:j+rw,k:k), fo, TCI=.false.)
                    q(pes(4):pee(4),i,j,k) = fo(rts:rte)
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            case (left)
              associate (q     => mesh%domain(ngb%ptr%id)%Jab  , &
                         recon => recon_arrays(ngb%ptr%id)%recon_hor)
              !$OMP PARALLEL DO PRIVATE(fo, ierr) COLLAPSE(3)
              do k = kds, kde
                do j = jts, jte
                  do i = its, ite
                    call recon(i,j)%ptr%reconstruct(q(pc,i-rw:i+rw,j-rw:j+rw,k:k), fo, TCI=.false.)
                    q(pes(1):pee(1),i,j,k) = fo(rls:rle)
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            end select
          end if
          end associate
        end do
      end if
    end do
    end associate

  end subroutine reconstruct_metric

  subroutine hope_recon_final()

    integer idom

    if (allocated(recon_3d_obj  )) deallocate(recon_3d_obj  )
    if (allocated(recon_hor_obj )) deallocate(recon_hor_obj )
    if (allocated(recon_ver_obj )) deallocate(recon_ver_obj )
    if (allocated(recon_topo_obj)) deallocate(recon_topo_obj)
    do idom = 1, mesh%ndom
      if( allocated(recon_arrays) )then
        if (allocated(recon_arrays(idom)%recon_3d  )) deallocate(recon_arrays(idom)%recon_3d  )
        if (allocated(recon_arrays(idom)%recon_hor )) deallocate(recon_arrays(idom)%recon_hor )
        if (allocated(recon_arrays(idom)%recon_ver )) deallocate(recon_arrays(idom)%recon_ver )
        if (allocated(recon_arrays(idom)%recon_topo)) deallocate(recon_arrays(idom)%recon_topo)
      endif
    end do

  end subroutine hope_recon_final
  
  function deconvolution(q,k,kds,kde)
   real(r_kind) :: deconvolution
   real(r_kind) :: q(5)
   integer :: k, kds, kde
   
   real(r_kind), dimension(3), parameter :: coef_mtxL = [ 23., 2., -1. ] / 24.
   real(r_kind), dimension(3), parameter :: coef_mtxC = [ -1., 26., -1. ] / 24.
   real(r_kind), dimension(3), parameter :: coef_mtxR = [ -1., 2., 23. ] / 24.
   
   if( k == kds )then
     deconvolution = dot_product( coef_mtxL, q(3:5) )
   elseif (k == kde )then
     deconvolution = dot_product( coef_mtxR, q(1:3) )
   else
     deconvolution = dot_product( coef_mtxC, q(2:4) )
   end if
  
  end function deconvolution
  
  !function convolution(q,k,kds,kde)
  ! real(r_kind) :: convolution
  ! real(r_kind) :: q(5)
  ! integer :: k, kds, kde
  ! 
  ! real(r_kind), dimension(3), parameter :: coef_mtxL = [ 25., -2., 1. ] / 24.
  ! real(r_kind), dimension(3), parameter :: coef_mtxC = [ 1., 22., 1. ] / 24.
  ! real(r_kind), dimension(3), parameter :: coef_mtxR = [ 1., -2., 25. ] / 24.
  ! 
  ! if( k == kds )then
  !  convolution = dot_product( coef_mtxL, q(3:5) )
  ! elseif (k == kde )then
  !   convolution = dot_product( coef_mtxR, q(1:3) )
  ! else
  !   convolution = dot_product( coef_mtxC, q(2:4) )
  ! end if
  !
  !end function convolution
  
  function convolution(q,k,kds,kde)
   real(r_kind) :: convolution
   real(r_kind) :: q(5)
   integer :: k, kds, kde
   
   real(r_kind), dimension(3), parameter :: coef_mtxL = [ 25., -2., 1. ] / 24.
   real(r_kind), dimension(5), parameter :: coef_mtxC = [ -17., 308., 5178., 308., -17. ] /5760.
   real(r_kind), dimension(3), parameter :: coef_mtxR = [ 1., -2., 25. ] / 24.
   
   if( k <= kds+1 )then
    convolution = dot_product( coef_mtxL, q(3:5) )
   elseif (k >= kde-1 )then
     convolution = dot_product( coef_mtxR, q(1:3) )
   else
     convolution = dot_product( coef_mtxC, q(1:5) )
   end if
  
  end function convolution

end module hope_recon_mod
